Discontinuous Galerkin methods for general-relativistic hydrodynamics: formulation and 
application to spherically symmetric spacetimes 



David Radice 

Max-Planck-Institut fur Gravitationsphysik, Albert Einstein Institut, Potsdam, Germany 

Luciano Rezzolla 

Max-Planck-Institut filr Gravitationsphysik, Albert Einstein Institut, Potsdam, Germany and 
Department of Physics and Astronomy, Louisiana State University, Baton Rouge, USA 

We have developed the formalism necessary to employ the discontinuous-Galerkin approach in general- 
relativistic hydrodynamics. The formalism is firstly presented in a general 4-dimensional setting and then 
specialized to the case of spherical symmetry within a 3 + 1 splitting of spacetime. As a direct application, 
we have constructed a one-dimensional code, EDGES, which has been used to asses the viability of these meth- 
ods via a series of tests involving highly relativistic flows in strong gravity. Our results show that discontinuous 
Galerkin methods are able not only to handle strong relativistic shock waves but, at the same time, to attain very 
high orders of accuracy and exponential convergence rates in smooth regions of the flow. Given these promising 
prospects and their affinity with a pseudospectral solution of the Einstein equations, discontinuous Galerkin 
methods could represent a new paradigm for the accurate numerical modelling in relativistic astrophysics. 

PACS numbers: 04.25.Dm, 02.70.Dh, 47. 11. Kb, 04.40.Dg, 



I. INTRODUCTION 

Special and general-relativistic hydrodynamics play a fun- 
damental role in a number of astrophysical scenarios char- 
acterized by strong gravitational fields and flows with high 
Lorentz factors. These scenarios are of central impor- 
tance for many phenomena in high-energy astrophysics or 
gravitational-wave astronomy, such as, core-collapse super- 
novae, coalescing neutron stars, accretion flows and relativis- 
tic jets. 

It is not surprising therefore that the study of numerical 
methods for relativistic hydrodynamics in astrophysical con- 
texts was an active branch of research already in the late '60s, 
starting with the seminal works by May and White [1] and 
by Wilson [2] and with the availability of the first massive 
computational facilities. In these investigations, the approach 
was to cast the relativistic hydrodynamics equations as non- 
linear advection equations in a form that resembles the New- 
tonian Euler equations. These were then solved using finite- 
difference (FD) schemes, stabilized using a combination of 
upwinding and artificial-viscosity methods to avoid excessive 
oscillations at shocks (see [3] for a comprehensive list of ref- 
erences). Although these methods allowed to perform the first 
numerical studies in general-relativistic hydrodynamics, they 
also had several limitations, such as the difficulty of tuning the 
artificial viscosity to avoid excessive smearing of the shock 
fronts or, most importantly, the limitation to mildly relativis- 
tic flows, i.e. with Lorentz factor W < 2 [3]. 

A major leap forward in numerical relativistic hydrodynam- 
ics took place when it was realized that the major problem be- 
hind Wilson's approach was the use of a formulation which 
breaks the conservative nature of the equations [4]. This real- 
ization lead to a formulation of the equations of relativistic hy- 
drodynamics in a conservative form, the so-called "Valencia 
formulation" [5], and to the use of finite-volume (FV) and FD 
high-resolution shock capturing (HRSC) methods for their nu- 



merical solution (see, e.g. [6, 7]). These methods were shown 
to be able to handle ultra-relativistic flows and to sharply re- 
solve shocks without spurious oscillations or need for artificial 
viscosity. For these reasons they have been the key ingredi- 
ent in a number of recent achievements of numerical relativis- 
tic hydrodynamics and magnetohydrodynamics (MHD; see, 
e.g. [8, 9] and references therein). 

Even tough FV and FD schemes have been particularly suc- 
cessful and are indeed the standard choice for modern numer- 
ical codes in relativistic hydrodynamics and MHD, they also 
suffer from some limitations, such as the difficulty of han- 
dling complicated grid structures and boundary conditions, 
or those associated with achieving high orders of accuracy. 
These are mainly due to the fact that high-order accuracy is 
generally attained with the use of large reconstruction stencils 
and expensive non-linear limiting operators, which become 
quickly cumbersome to handle when the grid is not struc- 
tured and/or quadratures are required in the computation of the 
fluxes, i.e. for higher than third order FV schemes. Large re- 
construction stencils also come with large ghost regions when 
doing parallel calculations, leading to poor scalability results. 
Finally these schemes are also often overly dissipative in situ- 
ations in which shock waves are not the dominant part of the 
dynamics and may fail to properly resolve fine structures of 
the flow [10]. This has important consequences for the accu- 
racy of general-relativistic hydrodynamics codes [11]. 

For these reasons, alternative approaches to general rela- 
tivistic hydrodynamics such as finite-element methods [12], 
or spectral methods [13, 14] are worth consideration. This lat- 
ter approach is particularly interesting because spectral meth- 
ods are able to attain very high accuracy, but it is also limited 
by the well known fact that these methods fail spectacularly 
when the solution develops large gradients or discontinuities. 
For this reason, spectral methods for relativistic hydrodynam- 
ics have been limited to the generation of initial data [15, 16] 
or to situations in which strong gradients could be treated 
with shock-tracking techniques within a multi-domain frame- 
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work [13]. 

More recently, however, a novel method was suggested 
by Dumbser and Zanotti [17], who presented a hybrid 
FV/discontinuous Galerkin (DG) approach for special rela- 
tivistic resistive magnetohydrodynamics (MHD). In this ap- 
proach a local spacetime DG method was used as an implicit 
predictor step in the context of a high-order FV scheme, in 
order to treat the stiff source term of resistive MHD. 

Indeed, Spectral Discontinuous Galerkin Methods (SDGM) 
and their variation employing Gaussian numerical integration 
(SDGM-NI) were developed to overcome some of the above 
limitations of FV and spectral or pseudospectral methods re- 
spectively [18]. These methods work essentially by com- 
bining the classical Runge-Kutta discontinuous Galerkin ap- 
proach by Cockburn [19] with the spectral element method 
(SEM) of Patera [20]. For this reason they are also often re- 
ferred to as DG-SEM or DG-SEM-NI. These methods are par- 
ticularly well suited for the solution of conservation laws and 
have been successfully applied to a number of classical hy- 
perbolic, parabolic and elliptic problems (see, e.g. [18, 21]). 
Finally they have been also successfully applied to the solu- 
tion of the Einstein equation in vacuum by Zumbush [22] and 
Field et al. [23]. 

We develop here the necessary formalism for the appli- 
cation of fully explicit DG methods to relativistic hydrody- 
namics on curved spacetimes. As an application we present 
a prototype code employing SDGM-NI for general relativis- 
tic hydrodynamics in spherical symmetry. We show that the 
proposed scheme is able to properly resolve strong shocks 
and achieve high-order, spectral accuracy for smooth solu- 
tions. While we will not discuss explicitly the coupling of 
the solution of the hydrodynamics equations with that of the 
Einstein equations, it is clear that a natural choice would be 
to use discontinuous Galerkin methods, such as the ones re- 
cently proposed by [22] or [23], or finite-element-methods 
such as the ones introduced by [24, 25], also for the metric 
evolution equations. This approach would have the advan- 
tage, with respect to the solution proposed with the "Mariage 
des Maillages" [26, 27], that the fluid and the spacetime vari- 
ables would share the same grid and no expensive interpola- 
tions would therefore be needed. 

The paper is organised as follows. In Sect. II we derive the 
general theory for the application of discontinuous Galerkin 
methods to relativistic hydrodynamics in curved spacetimes 
and we specialize it to the spherically-symmetric case. In 
Sect. Ill we present our prototype numerical code, EDGES 
(Extensible Discontinuous GalErkin Spectral library), which 
was used to test DG methods for general relativistic hydrody- 
namics in one-dimension (ID) and spherical symmetry. The 
results obtained on a representative number of test cases are 
then presented in Sect. IV. Finally Sect. V is dedicated to the 
summary and conclusions. 

We use a spacetime signature (— , +, +, +), with Greek in- 
dices running from to 3 and the Latin indices from 1 to 3. 
We also employ the standard convention for the summation 
over repeated indices. Finally, all the quantities are expressed 
in a system of units in which c = G = M Q = 1, unless 
otherwise stated. 



II. DISCONTINUOUS GALERKIN METHODS FOR 
GENERAL-RELATIVISTIC HYDRODYNAMICS 



Broadly speaking, Galerkin methods are projection meth- 
ods for the weak formulation of the equations. In the case 
of the general-relativistic hydrodynamics equations, such for- 
mulation could be obtained in two different ways. The first 
one consists in starting from the relativistic hydrodynamics 
equations written in a conservative form in a chosen coordi- 
nate system, e.g. the Valencia formulation [5], and then in- 
tegrating them against a test function. The numerical scheme 
obtained with the Galerkin projection would then be a direct 
generalization of the standard HRSC schemes used in general- 
relativistic hydrodynamics. Indeed, when considered at first- 
order only, DG schemes reduce to FV ones and it is for this 
reason that they are often interpreted as an alternative way to 
attain high-order FV methods 



While this approach is certainly possible and would seem 
to be quite natural, it has the limitation that since we start 
from the equations in their coordinate form, we also have to 
choose a metric with respect to which the volume integrals are 
performed. The choice of such a metric is effectively arbitrary, 
but any choice different from that of a flat metric corresponds 
to the absorption of a multiplicative factor into the definition 
of the test function and thus it is equivalent to a modification 
of the scheme in the higher-than-first-order case. As a result, 
the choice of the metric is unimportant only in the FV limit of 
Galerkin methods, but it plays a central role in higher-order 
Galerkin schemes. 



A second way to obtain the weak formulation of the equa- 
tions and the one actually outlined in this paper, is to follow an 
approach which is instead manifestly covariant and thus does 
not require any assumed background. After the formulation 
is obtained, it can then be decomposed in the standard 3 + 1 
split of general relativity. This choice has the advantage of 
producing the most natural extension of the commonly used 
HRSC frameworks to the DG case. The resulting schemes 
will be naturally covariant, suited for the use with standard 
spacelike or null foliations or even independently of any fo- 
liation or coordinate system. The reason why this is possible 
lays, as pointed out by Meier [12], in the covariant nature of 
finite-element methods and, by extension, of DG methods. In 
these methods, in fact, the equations are formulated on refer- 
ence elements mapped into the physical space via diffeomor- 
phisms, thus removing any need for a (preferred) coordinate 
system. The important difference between our approach and 
the one by Meier [12] is in the use of non-conforming, dis- 
continuous, Galerkin methods. This gives us the possibility 
of reducing the coupling of the numerical solution across the 
elements to flux terms, thus enabling the construction of glob- 
ally explicit, local schemes, in contrast to the need of solving 
implicit, global, nonlinear problems. 
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A. Weak formulation of the equations of relativistic 
hydrodynamics 

Let (Ai,g a p) be a strongly hyperbolic, C 2 , spacetime with 
metric g a p and let V be the covariant derivative associated 
with g a p. We consider a perfect fluid described by a rest- 
mass-density 4-vector J a and a stress energy tensor T Q/3 de- 
fined by 



J a = pu a , 



phu a vr 



pg 



a/3 



(1) 



where p is the rest-mass density, u a is the fluid 4-velocity, 
p is the pressure, e is the specific internal energy and h = 
1 + e + p/pis the specific enthalpy. 

If we assume baryon-number conservation and a generic 
equation of state (EOS) of the form p — p(p, e), then the 
equations of motion for the fluid on M. read 

v a j a = o, v^ = o, p = p (p,s). (2) 

In general these equations are to be intended in the sense of 
distributions 1 , since we expect the solution to develop singu- 
larities in the form of shock waves. 

In general we are interested in solving (2) on an open, reg- 
ular 2 , finite domain O C M, with suitable initial/boundary 
data. A precise mathematical formulation of this problem can 
be done within the context of bounded divergence-measure 
vector fields using the theory developed in [28, 29]. In partic- 
ular, we will look for solutions in the functional space, V, of 
all the i.e. "bounded", vector fields over il, whose diver- 
gence, in the sense of the distributions, are Radon measures 3 

As a first step we introduce a triangulation of N "elements " 
of ft, {Qj}f=i, by selecting a family of diffeomorphisms 
tfj : K C R 4 -» O, Slj = (fij (K) such that 



N 

3=1 



n Q.j = 



(3) 



where K is the, so-called, "reference element", usually an 
hypercube or a 4D simplex and Clj denotes the interior of flj. 
We also arrange the local coordinate system, {x^}j, induced 
by <pj, so that 9 is timelike or null. 

If we now look for solutions J a e V, the first of the equa- 
tions (2) is equivalent, in the sense of distributions, to 



JV 

£ 

3=1 



[ J a \7 a <t>e- [ <j)J a e aM5 ^Ax^Ax s 



0. 



(4) 



1 Here and hereafter we will denominate as distributions generalized func- 
tions, such as Dirac delta. 

2 See [28] for a detailed discussion of the regularity requirements. Broadly 
speaking this amounts to having a domain which has a normal defined ev- 
erywhere except for at most a discrete set of points (vertexes); i.e. a cubic 
box is a regular domain. 

3 For a precise definition see [29]. Broadly speaking this condition means 
that we restrict ourselves to cases in which the solution presents at most 
mild singularities, such as jump-discontinuities. 



all (j) e CUM) 4 . Note that in 
symbol e refers to the proper 
the spacetime, i.e. in any 



local 



expression (4) 
volume form 
chart, {x^}, 



for 
the 
of 

e = y^ge A e 1 A e 2 A e 3 = e a p 7 s<ix a dx l3 dx~ f dx 5 , and 
J a is the internal normal trace of J a . This object reduces 
simply to J a , when J a and ilk are regular, but in the general 
case the second integral has to be intended as the action of a 
measure, J a € a/ 3 1 sdx l3 dx' < dx s , on <j> [28]. 

In the same way, if we look for solutions T a/3 e V ® V, 
the second of the equations (2) is equivalent, in the sense of 
distributions, to 



N 

„ , JdQi 



(5) 



3~1 3 



for all the one- forms <p a e Cl(M;T*M), T*M being the 
co-tangent bundle of A4. Again, T a ^ is a generalization of 
T af) and the integral has to be interpreted as the action of 
T a Pe l 3 1 s t j,dx~ < dx s dx fJ - on <f) a in the non-smooth case. 

The solution of the relativistic hydrodynamics equations 
consists then in finding 

feV s.t. (4) holds V^eCo^M), (6a) 

T afi g V®V s.t. (5) holds \/cj) a g Cl(M;T*M), (6b) 

together with an EOS and proper boundary-initial data, to be 
specified through J a and T a ^ on dCl. We remark, again, 
that (6) is perfectly equivalent, in the sense of distributions, 
to (2) and that the triangulation has been introduced 

mainly for later convenience. 



B. Spacetime discontinuous Galerkin formulation 

As mentioned above, within the Galerkin approach a nu- 
merical scheme is obtained by projecting (6) on a finite di- 
mensional subspaces V C V. In general this space is con- 
structed starting from the space of piecewise polynomials, in 
particular we define 

X = {u G L°°(n) : u o ipj G V D (K), j = l,...,N}, (7) 

where ¥]j(K) is the space of polynomials with at most degree 
D on K. Notice that the functions in X are allowed to be dis- 
continuous at the edges of the elements, hence the name "dis- 
continuous Galerkin" for the resulting numerical scheme. The 
space V is taken as the space of all the vector fields "whose 
components are elements of X", more precisely 

V={u a eV: [^}*u a e [Fd(K)]\ j = 1,...,N}, (8) 



4 We recall that a function of class Cq (H) is a function of class C" (fi) and, 
in addition, with compact support in CI. 



4 



where pP D (Jf)f is the space of 4-tuples of polynomials with 
at most degree D on K and [<fj]* is the pull-back associated 
with <pj. The Galerkin method is then simply the restric- 
tion of (6) to V so that it consists in finding J a g V and 
T afi eV ®V such that (4) and (5) hold for suitable choices 
of the test functions, <fi and (f> a , that we are going to discuss in 
the following. The resulting equations can be solved numer- 
ically, because a finite number of conditions suffice to fully 
determine J a and T a @ as long as we have a way to evaluate 
the fluxes J a and T a/3 on dSlj . 

As discussed above, in the continuous case these fluxes are 
simply the restriction of J a and T a/3 to dSlj . More explicitly 
if we write symbolically V = {p, v},u 2 ,u 3 , e} for the prim- 
itives variables and consider J a and T Q/3 as functions of V, 
i.e. J a = J a {V) and T a P = T a ^(V), then J a = J a (V*) 
and T af} = T afj (V*), V* being the restriction of V on dSlj, 
as the fluxes can only depend on the location in the space- 
time through V. Stated differently, J a and T a/3 are simply 
the Godunov fluxes for the conservation law. In the general 
case J a and T a ^ can be determined with causality consid- 
erations on spacelike boundaries 5 , or as solutions of general- 
ized Riemann problems, on timelike and null-like boundaries, 
as they are known as soon as V* is known on those bound- 
aries. This is basically analogous to the Newtonian case when 
using spacetime DG methods to discretize balance laws. For 
this reason we refer to [30] for a more in-depth discussion. 
We limit ourselves to note that, in the context of a numerical 
scheme, the computation of the fluxes can be greatly simpli- 
fied with the use of approximate Riemann solvers, such as the 
HLLE scheme. In that case J a and T a ^ are approximated 
directly without the need to explicitly compute the solution of 
the Riemann problem, V*, at the interface. This is discussed 
in some more detail, for the 3+1 case, in the Appendix B. 

Once we have a way to compute the fluxes, the fully dis- 
crete equations are readily obtained by testing (4) and (5) on 
a set of linearly-independent test functions, <f> and <f> a . The 
resulting finite set of equations can be cast in a set of non- 
linear equations for the spectral coefficients of the numerical 
solutions, when these are expanded over a linear basis of V, 
or V <8> V, for J a and T al3 respectively, and by following a 
standard finite-element method procedure, see e.g. [31]. In 
particular, it is sufficient to consider test functions <j> g X and 
4> a g V, so that the final numerical scheme consists in finding 



r g v 

T af) g V <8> V 



s.t. (4) holds 
s.t. (5) holds 



el, 
a g V . 



(9a) 
(9b) 



where the test functions on dflj in the boundary integrals ap- 
pearing in (9) has to be interpreted as being a Cq extension to 
M. of the original test function, created in such a way as to 



smoothly match the one-sided limit, from the interior of Clj, 
of the original test function. 

To explicitly write down the method, in every finite el- 
ement, ilj, we select a set of conserved quantities, C = 
{J°,T°' i } 6 , for which there exists a one-to-one relation, in- 
volving the EOS, with the set of primitive variables [32], V, 
so that we can formally write J 1 = J\C) and = T^(C). 
We can then obtain a set of nonlinear equations for C g X 5 
simply expanding the Galerkin conditions (4) and (5) 



N 

E 

3 = 1 



[ j°a o 0e+ / j\c)d 3 ^ 

N 

=1 Jdn 3 



(10) 



and, setting <f> a — (j)d' 1 a , 



N 



E 



/ T°»d cbe+ f T^iQdrf 

f T^0e 1 , Q/37 dx Q dx /3 dx''+ (11) 
Jon, 



N 

E 

3 = 1 



where T a j3l are the Christoffel symbols and <f> g X. 

The key point here is that, as the functions are discontinu- 
ous across the dClj's, these equations are local equations for 
the spectral coefficients within the Oj 's coupled only through 
the fluxes. In particular this implies that, if the computa- 
tional grid "follows the causal structure of the spacetime", 
in the sense that it can be traversed with a succession of 
"causal slices" satisfying some sort of generalized Courant- 
Friedrichs-Lewy (CFL) condition ensuring the causal discon- 
nection between the timelike boundaries of the elements, then 
the discontinuous Galerkin method becomes globally explicit. 
Under these conditions, in fact, the fluxes between the ele- 
ments of the grid slice depends only on the data on the pre- 
vious slice. Once these are computed, we are left with a set 
of formally decoupled equations involving the spectral coef- 
ficients of the numerical solution in the different elements. 
The precise mathematical definitions of "causal slices" and 
of "grid that follows the causal structure of the spacetime" are 
given in Appendix A. 



C. Discontinuous Galerkin formulation in the 3 + 1 split 

Clearly the strategy outlined above could be useful in situa- 
tions in which stiff sources are present, such as in the resistive 



5 Where the solution has different limits, V\ and V2 from different sides of 
T C dtlj and T is spacelike, we proceed as in the Godunov method and 
we set J a = J a [P* (Pi ,V 2 )] - Causality requires that V* must depend 
only on the past limit of V at T, say Vi, this implies V* = Pi, thus 
J" = [J a ]i - The same argument can be applied to evaluate T al3 . 



6 Other choices are possible, for example in the context of a 3 + 1 split we 
could use the same conserved quantities as the ones used in the Valencia 
formulation. 
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MHD case, where implicit DG methods have already been 
shown to be suitable as predictors to treat stiff sources [17], 
but in the unmagnetized case an implicit time stepping is un- 
necessarily expensive in most situations. The generation of a 
triangulation which follows the causal structure of the space- 
time could also be highly non-trivial, especially if the space- 
time is evolved dynamically. For these reasons, instead of 
directly solving (10) and (11), one can use them in order to 
derive a fully explicit scheme. This can be accomplished by 
performing a 3 + 1 split directly at the level of the discrete 
problem (9). 

As customary, we foliate the spacetime along t = const, hy- 
persurfaces, E t , and consider a vector t a such that t a V Q i = 
1. Using this vector we define the three-volume form rj a p 7 — 
t&ap-yt 5 ■ Also, as usually done in this context, we can use 
the integral lines of t a to identify points on E t with points on 
E = Eo and interpret the variation of the fields across the E t 's 
as being the result of a dynamics on a three-manifold, E. We 
are then interested in studying (2) in a world-tube S x (0, t), 
ScS being an open, bounded, regular domain in E, together 
with proper boundary-initial conditions. Clearly this is a par- 
ticular case of the general problem studied above. 

In order to apply the discontinuous Galerkin formulation, 
we consider a triangulation {Sj}^^ of S, by selecting a fam- 



ily of diffeomorphisms ■ : T C 
that 



E, Sj = $j(T) such 



N 



3 = 1 



Si n 



(12) 



where T is now a three-dimensional reference element, a cube 
or a tetrahedron. This induces a triangulation of £1 

at n < ^ N 'Q 

{^>}f=?,„=i = { S i x ( nAt > ( n + !) A *) . , t > ( 13 ) 

that follows the causal structure of the spacetime on the n-th 
thin-sandwich fl n = S x (t n , t n + At), at least for small At. 

As we intend to use the method of lines in order to integrate 
the equations in time, we can factor out the time dependence 
from the test functions by considering the functional space 

Y = {u G L°°(S) : uo*,- e ¥ D (T), j = l,...,N}, (14) 

and vector functions "whose components are elements of Y": 

W={u a eV: l$jlu a e [F D (T)]\ j = l,...,N}. 

(15) 

Given a function iieF.we can consider it as a function, u, 
over S x (0, t) with the identification u(x l , s) = u(x l ) so that 
Y <-> Y C X and, with a slight abuse of the notation, we 
can consider Y as being a subspace of X. In a similar way 
we can consider W to be a subspace of V. For these reasons 
we can choose e Y in (10) and (1 1). We can thus obtain a 
fully explicit method by re-projecting (9) onto a new couple 
of subspaces or, equivalently, by finding 

J a G C 1 ((0, t); W) s.t. (4) holds G Y , (16a) 

T alB £ ^((O^W^W) s.t. (5) holds V^ a G W. (16b) 



In the expressions above, and as customary when dealing 
with evolution problems, we have used the notation u(-) G 
C fc ((o, b); X), where a, b G K and X is a Banach space, to 
indicate that, when u is a regarded as a function only of time, 
it describes a regular, C k , curve in X. In other words, when 
u(t, x\) is interpreted as a function of time only, it is of class 
C k , while when u(ti,x l ) is interpreted as a function of x % 
only, it is an element of the function space X. Note also that 
in this new formulation we do not allow the numerical solu- 
tion to be discontinuous in time, so that the time integration 
can be performed with a standard solver for ordinary differen- 
tial equations (ODEs). 

We can derive a more explicit form for the (16) by project- 
ing (10) with (f> G Y, dividing both terms by At and by letting 
At -> 0, to obtain 



jv . jv r . 

j = l j = l \_-' S 3 



/ T 4>r) ia pdx a dx li 

JdS, 



(17) 



Reasoning along the same lines, we can derive an explicit 
discretization of the second of the (2), starting from the (11), 
to obtain 



j = l , ' S i j=l L Sj 

[ T l ^m^dx a dx^ - [ T a P(C)r» a <f>ri 



(18) 



Finally the 3 + 1 discontinuous Galerkin formulation can be 
summarized as in finding 

J a G ^((O,*); W) s.t. (17) holds V0GF, (19a) 

T af} G ^((O.t); W®W) s.t. (18) holds V0" G W. (19b) 

This scheme can be interpreted as a higher-order general- 
ization of a FV discretization of the manifestly covariant for- 
mulation of relativistic hydrodynamics proposed by [32]. As 
a consequence, our scheme inherits properties such as hyper- 
bolicity and the flexibility to work with spacelike or null-like 
foliations directly from [32]. This can be seen considering the 
case in which D — 0, that is, when looking for solutions that 
are constant over each element, Sj . Then a sufficient number 
of Galerkin conditions can be obtained by simply choosing 
(f> = xsj for j = 1, 2, . . . , N, where \e is the indicator func- 
tion of the set E, i.e. a function which is equal to one in E and 
identically zero elsewhere. With this choice we obtain the set 
of equations 

dt I J°V+ f J l (C) Via pdx a dxP =0, (20) 

JSi JdS, 



and 



■ / T^7)+ f 

J Si Jdi 



d t I T^r,+ T^^ afS dx a dx p = 

IdSi 



J S-j 



(21) 
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for all j = 1,2, ... ,N. These can easily be recognised as 
being the FV discretization of the formulation by [32]. 



D. Discontinuous Galerkin formulation in spherical symmetry 

As a particular case of the formalism outlined above we 
consider the case in which the spacetime is spherically sym- 
metric. This case is particularly interesting because the equa- 
tions become 1 + 1 dimensional and are therefore well-suited 
for rapid prototyping and testing of new methods and tech- 
niques 7 . 

In particular we consider a spherically symmetric space- 
time in a radial-polar gauge 

,2 i,2 , a2 a 2 



factor. Furthermore we introduce the "conserved" quantities 

V = aAJ 1 = pAW , (24a) 

S = aT tr = phW 2 v , (24b) 

E = a 2 T u = phW 2 - p , (24c) 

t = E — V . (24d) 

With these definitions, the Einstein equations reduce to the 
Hamiltonian constraint 



d r m = Airr 2 E , 
and the slicing condition d t K eg = K de = 



d r v = A 2 



+ 4irr(p + Sv) 



(25) 



(26) 



ds = -a dt + A dr + r (d9 + sin 9d<j) ) , (22) where K is the extrinsic curvature. 



where a and A are functions of t and r only. We next intro- 
duce the Bondi mass function, m, and the metric potential, v, 
by 



i- 2m(t ' r) V V2 , a(t,r) = e^K (23) 



A(t,r) 



Following [37], we define the physical velocity, w, by jj = 
Au r /au f , where W — au l = (1 — v 2 )^ 1 / 2 is the Lorentz 



These equations have to be integrated with the boundary 
conditions given by m(0) = and by the requirement that v 
matches the Schwarzschild solution at the outer boundary of 
the computational domain (see e.g. [38] for a detailed deriva- 
tion and discussion of this equation). 

The equations for the hydrodynamics are simply (17) and 
(18), where the elements, Sj, are taken to be the spherical 



shells Ti < r < r 



■j+v 



These can be written in terms of the 



conserved quantities (24) and specialized for the metric (22) 
by substituting the explicit expression for the Christoffel sym- 
bols and the determinant of the metric, to obtain 



± r 

.1=1 Jr i 



N 



d t F t (C)<f ) r 2 dr^Y,\ f 



r 3 + l 



XF r (C)d r (f>r 2 dr - 



r 2 XF r , 



rj+i P+ 1 
r + 



s(C)(f>r 2 dr 



(27) 



where we define X = a /A. Here, the fluxes are given by 

F* = {V, S, t} , F r = {Vv, Sv+p, S- Vv} , (28) 
while the source term is 

s = I 0, {Sv-t-V) (saAitrp + aA ™ ) 

1 ' . (29) 

. m „ap n \ 

In the derivation of (27) the momentum constraint was used to 
substitute the derivatives of the metric in the source term and a 



Of course in 1 + 1 dimensions a Lagrangian approach such as the one 
presented in [33, 34] is by far superior as it allows for natural spatial adap- 
tivity and conservation properties. However, both of these advantages dis- 
appear in more than one spatial dimension. On the other hand Lagrangian 
approaches to multidimensional relativistic hydrodynamics have been re- 
cently proposed within the context of smoothed-particle-hydrodynamics 
schemes in special [35] and general relativity [36]. 



factor X was absorbed into the test function in the derivation 
of the equation for r. 

A close examination of the (27) reveals that, again, this for- 
mulation of the equations can be interpreted as an higher-order 
generalization of a classical FV discretization of the equations 
of relativistic hydrodynamics. In particular it can be seen that, 
in the case in which V , S and r are constant over each ele- 
ment, (27) reduces to the FV method discussed in [37]. 



III. THE EDGES CODE 

In order to test discontinuous Galerkin methods for rela- 
tivistic hydrodynamics and their reduction to spherical sym- 
metry, we have developed a new ID code, EDGES. This con- 
sists of a general discontinuous Galerkin library which is then 
used by a code solving the general-relativistic hydrodynam- 
ics equations in spherical symmetry. EDGES makes exten- 
sive use of advanced generic programming techniques such as 
static polymorphism via recursive templates, and expression 
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templates see e.g. [39]. The code employs the Blitz++ 
high-performance array library [40, 41] and makes use of the 
UMFPACK multi-frontal sparse factorization method [42-46] 
for linear systems inversion. In what follows we describe in 
more detail the implementation of the discontinuous Galerkin 
approach within EDGES. 



A. The DG equations in a fully discrete form 

We consider a spherical shell S — [0, R] in the spacetime 
(22), containing a fluid described by (1) and (2). Furthermore 
we consider a "triangulation" of S = UjLi &j wnere ( see 
Fig. 1 for a scheme of the triangulation) 



Si n Sj = , Vi/j, Sj = <Pj([-l, 1]) • (30) 
The functional space that we consider in EDGES is, 

Z={nei M (S):«o W eP l) [-i,i]}, (31) 

the set of all the functions that are polynomials of degree D 
over each Sj . 

If we denote by Lu(x) the D-th Legendre polynomial on 
[—1, 1], a Gaussian quadrature of order 2D— I can be obtained 
with the formula 

-1 D 

/ /(iJdiw^Wj/W, (32) 

J- 1 i=0 

where {xi}fL are the zeros of (1 — x 2 )dL D (x) /dx, Wi are a 
set of weights given by 



w. 



= j li(x)dx, i = Q,l,...,D 



(33) 



and {k(x)}fL are the Lagrange polynomials associated with 
the nodes {xi}fL , i.e. a set of polynomials of degree D such 
that k(xk) — 5ik for i, k = 0, 1, ... , D. Given two regular 
functions / and g in r € (0, R) we now define their continu- 
ous scalar product as 

(f,g)= / f(r)g(r)r 2 dr, (34) 
Jo 

and use the quadrature formula (32) to introduce the discrete 
scalar product 



N 



(f,g)o = J2(f,g)j,D ■ 

3=1 



(35) 



where 



D 



if,9)j,D = w * Iv'jl f[<Pj( x i)] 9[<Pj{xi)] (36) 



i=0 



and \<p'j\ is the Jacobian of the affine transformation 

<Pj- [-1,1] -> Sj. 




Figure 1. Scheme of the spacetime grid structure in EDGES. The col- 
location nodes (filled squares) are generated and stored on the refer- 
ence element, T, and then mapped onto the finite element Sj with an 
affine transformation, ipj. Note that each point on the boundary of 
an element is associated with two distinct degrees of freedom, thus 
allowing the functions to have two different one-sided limits. Be- 
cause the number of collocation points needed is D + 1, where D 
is the order of the polynomial representation, the figure refers to a 
polynomial of order three. 



With these definitions in place we can construct a fully dis- 
crete systems by looking for solutions T>, S,t € Z and com- 
puting the integrals in (27) using the Gaussian quadrature (32) 
over each element. In particular, using the notation (35), we 
obtain 



(r 2 d t F\cf) D = (r 2 XF r ,d r < 

N 



n 



r 2 XT r c 



Id 

rj+i 



+ ( 



r s, 



(37) 



' D ' 



where rj is the left vertex of the j-th element in the discretiza- 
tion of the radial coordinate in the line element (22). 

Considering now test functions with support contained 
within a given element, Sj, say 



cj>(r) = U(<p. 1 {r))x[r ] ,r ]+1 \{r) 



(38) 



where k (ip- 1 (r)) is the Lagrange polynomial or order i eval- 
uated at the position x = ipj 1 ^) of the reference element 
which is mapped into r, while X[r j ,r j+1 ] is the indicator func- 
tion and thus equal to one in the interval [rj,r J+1 ] and zero 
elsewhere. Expanding F* over the Lagrange basis of Sj 

D 

[F t o i p j ](x) = J2F t jk l k (x), (39) 

fc=0 

we obtain a set of coupled ordinary differential equations for 
the coefficients F\ ,• 



r%d t F% = (r 2 XFr, dr k) jD 
1 



r 2 XT r h 



r 3 + l 



(40) 



where = (fij(xi) and we used the fact that (k,lk)j,D = 
Sik- An explicit example of (40) for the continuity equation 
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is given in the Appendix B. Hereafter we will consider only 
affine maps ipj, so that the transformation between the refer- 
ence element T and the finite element S is given by 



iPj(x) 



1 



rj+i 



1 



(41) 



and |vJj-| = \r j+1 - rj\/2. 

Figure 1 offers a schematic representation of the spacetime 
grid structure in EDGES. The D + 1 collocation nodes (i.e. the 
four filled squares for a polynomial of degree three) are gener- 
ated and stored on the reference element, T, and then mapped 
onto the finite element Sj with an affine transformation, (pj. 
Each point on the boundary of an element is associated with 
two distinct degrees of freedom, so that the corresponding val- 
ues of the functions can be different there. 

This system of equations is coupled with the equations for 
the evolution of the metric quantities, described in Sect. Ill B, 
and closed with an EOS. In particular we use an ideal-fluid (or 
T-law) EOS 



p=(T-l)pe. 



(42) 



Expression (40) is clearly a collocation scheme for V, S 
and r on the grid illustrated in Fig. 1, since we can inter- 
pret the expansion coefficients as collocated values, i.e. fji = 
f(rji) for any function /. In our code we evaluate the fluxes, 
T r , using the relativistic HLLE approximate Riemann solver 
and evolve numerically the resulting set of equations using a 
second-order strongly stability preserving (SSP) Runge-Kutta 
method [47]. We note that, as expected, if D = 0, then 



U 



const., the first term on the right-hand side vanishes 



and we are left with a standard FV scheme. 

A final note concerns the CFL condition needed to ensure 
the linear stability of the scheme. The linear stability of Leg- 
endre pseudospectral methods for hyperbolic equations has 
been studied in [48], where it has been shown that instability 
can be obtained if At ~ D~ 2 . As a result, in our code, we 
use a timestep given by 



At 



a 



CFL 



(D + iy 



(43) 



where Ax is the size of the smallest element, the Ccfl is a co- 
efficient that is reminiscent (but distinct from) the traditional 
"CFL factor" and has to be determined empirically. This co- 
efficient is usually taken to be Ccfl ~ 0.2 — 0.3, but our 
numerical experience (at least in ID) seems to suggest that 
this condition is overly restrictive as stable evolutions can be 
obtained with Ccfl ~ 1 in most situations. 



B. Coupling with the spacetime 

As discussed in Sect. II D, in spherical symmetry and with 
the gauge chosen, the Einstein equations are simple con- 
straints on each spatial slice and thus ODEs in the form 



d r u(t,r) = f(t,r) 



which could be easily integrated to very high accuracy. How- 
ever, instead of using a standard Runge-Kutta method for their 
integration, we found that a more efficient and accurate ap- 
proach to solve equations (25) and (26) consists in using an 
implicit discontinuous spectral ODE solver which makes use 
of the same grid as the hydrodynamical variables. Such ap- 
proach has the advantage of obtaining a numerical solution 
with a degree of accuracy which is of the same order or higher 
with respect to the one attained for the hydrodynamical vari- 
ables, without requiring a very small step, as would have been 
the case for a Runge-Kutta method. 

In order to implement this approach it is necessary to al- 
low the solution to be discontinuous across the elements, 
while imposing the continuity using an interior penalty tech- 
nique [49, 50]. In particular the discontinuous Galerkin for- 
mulation that we use reads 



N 

E 

.?=i 



d r u (fidr — Hj\u\ 



N 

E 



r r j+i 

J Ta 



where <f>j = (j>(rf), 



f<pdr, 
(45) 

u(t,Tj ) — u(t,r^) is the jump 



term, and fj,j <~ l/Ar, is the penalization coefficient. It is 
straightforward to see that this penalization term, which we 
refer to as "upwind penalization" in contrast to the usual sym- 
metric penalization term J2f=i MjMjMj> nas tne effect of 
enforcing u(t, r+) = u(t, rj) without affecting the equation 

for the collocation value in rj . 

As customary in the context of DG-SEM-NI methods, we 
approximate (45) as 



N 



(d r u, 4>) D - ^2 MjHj 4>j = if, 4>)d ■ 



(46) 



so that if we arrange the values of u(t, r) and f(t, r) on the 
collocation points in two arrays u(t) and f(t), the previous 
can be written as 



Au(t)=f(t) 



(47) 



(44) 



where A is a large-sparse matrix. As A does not depend 
on t, in EDGES this matrix is pre-computed, stored and pre- 
factorized using UMFPACK (see e.g. [51]). At each step, then, 
we have simply to compute f(t) and use the factorized version 
of A to efficiently compute u(t) = j4 _1 /(t). 



C. Limiters, spectral viscosity and spectral filtering 

It is well known that high-order numerical methods suf- 
fer from numerical oscillations in the presence of disconti- 
nuities (Gibbs phenomenon). If these oscillations are not suit- 
ably handled, they tend to grow out of control and destabi- 
lize the method. To overcome this difficulty several different 
"flattening techniques" have been developed in the context of 
FD, FV and spectral methods to artificially lower the order of 
the methods in the presence of shock waves. Some examples 
are artificial viscosity, flux limiting, PPM or ENOAVENO re- 
construction. Many of these techniques are implemented in 
EDGES and can be activated during the evolution. 
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1. Complications of spherical symmetry 

While discontinuous Galerkin methods can in general be 
used in combination with the large majority of flattening tech- 
niques, this is not the case in spherical symmetry. The reason 
for this is that if we consider a function u and interpret any of 
the flatting methods above as the application of an operator C 
to u, it can be shown that, for the above mentioned flattening 
techniques, 

/ udx = C(u)dx, (48) 

which is the desired behaviour for a conservative scheme in 
ID. However, in the case of spherical symmetry, the conser- 
vation property that we should satisfy reads 

/ ur 2 dr = C(u)r 2 dr. (49) 
Jnj Jrtj 

Unfortunately, the property (49) does not hold for almost 
all of the stabilization techniques which will be discussed 
in Sect. IIIC2 and which are generally coordinate depen- 
dent (the only exception being given by the spectral viscos- 
ity method which is inherently a differential operator). In the 
case of FV codes this is not a problem because the reconstruc- 
tion is only used in the computation of the fluxes and the vol- 
ume averages are evolved using the correct cell volume, but 
in a scheme that works with the actual point-wise value of the 
solution, this leads to unacceptable variations of the volume 
integral of u (e.g. of the total "mass") in the elements with 
r <C 1 and r»l, which in turn leads to the development of 
large numerical errors and/or instabilities. For this reason it is 
necessary to introduce a correction to the flattening procedure 
in the spherically symmetric case. 

A naive way to obtain the wanted result would be to modify 
the operator C as 

C(u) = ^C(ur 2 ) , (50) 

so that C now satisfies (49). In practice, however, this strategy 
results, as can be easily foreseen, in large numerical errors 
near r = and is this of little use. For this reason we had 
to adopt a more radical approach and add a "correction step" 
after the application of C to enforce (49). In particular we im- 
plemented three different strategies which we discuss below. 

The first one, which we refer to as "dummy" correction, 
consists in simply adding to C{u) a function C, constant on 
every element, such that (49) is satisfied. This is a very sim- 
ple approach and has the advantage of not increasing the total 
variation of C{u) inside the single elements, as conservation 
is basically obtained at the expense of the amplification of the 
jumps of C(u) across the elements. In this way the additional 
total variation generated by the corrective procedure is con- 
centrated over the grid points which constitute the "finite vol- 
ume part" of the method. 

The second one, which we refer to as "bubble" correction, 
consists in adding to C{u) a function b, which is a bubble 



function over each element, i.e. a function which is zero at the 
edges of the element, such as 

[bo l p j ](x)=K j (l-x 2 ), (51) 

where the Kj 's are chosen in each element so that (49) is satis- 
fied. This approach has the advantage of avoiding the creation 
of artificial discontinuities at the boundary of the elements, at 
the price of a small increase in the total variation of C(u). 

Finally the third one, which we refer to "intrinsic" correc- 
tion, consists in modifying the action of C so that L q = C(L q ) 
for q < 2. In this way both (48) and (49) are satisfied. In prac- 
tice this is obtained by overwriting the low-order coefficients 
of the discrete Legendre transform (DLT) of C(u) with the 
ones of u. This method has the advantage of not introducing 
any unwanted extra total variation and, for this reason, is the 
one that has the best mathematical properties. The only lim- 
itation is that it effectively weakens C and could thus fail to 
completely removing the Gibbs oscillations. 

All these three techniques can be basically used in combina- 
tion with any of the stabilization methods which are described 
in the following Section. 



2. Stabilization techniques 

The most commonly used stabilization technique in discon- 
tinuous Galerkin methods is based on slope limiting. These 
methods were originally developed for FV schemes [6], but 
are easily modified to work with discontinuous Galerkin 
schemes as done by Cockburn and Shu [19], who introduced 
the "Alii" limiter based on a generalized version of the "min- 
mod limiter". In our code we implemented a refined version 
of this limiter originally proposed by [52] and subsequently 
improved and extended by [53]. This essentially works by 
recursively limiting the coefficients of the spectral represen- 
tation of the solution on the various elements. The main ad- 
vantage of this technique is that it does not require tuning and 
it is usually very reliable, while its main limitation is that its 
use often results in excessive flattening of the solution in the 
presence of discontinuities. Moreover we found that all these 
methods perform rather poorly in conjunction with the cor- 
rection techniques outlined above and for this reason we have 
rarely used them. 

Another method implemented in EDGES is the "spectral 
viscosity" method first proposed by Maday et al. [54] in the 
context of Legendre pseudospectral methods for conservation 
laws. This method consists in the addiction to the right hand 
side of (40) of a term 

-e D (Q 3^,3^)0, (52) 

where Ed is a coefficient depending on the number of grid 
points and Q is a viscosity kernel whose action on a scalar 
function u reads, in every element, 

D 

[Qu](x) = ^2Q k u k L k (x) . (53) 

fc=0 



10 



Here, iik are the coefficients of the DLT of u and Qk are real 
numbers such that 



= 0. 



< Q k < 1 , 



for k < mo , (54a) 
for k > 772 £j . (54b) 



Note that mu effectively plays the role of the cut-off fre- 
quency of the filter. Meday et al. [54], were able to show that, 
in the context of scalar conservation laws in ID and using only 
one domain, this method is able to stabilize the scheme, with- 
out spoiling its spectral accuracy if 

e D ~ D~ e , for D > 1 , (55) 

where theta is just an exponent such that < 9 < 1 and 

m D ~ D q/4 , for < q < 9 < 1 . (56) 

In our code we have set 

At 1 



e D = /i 



N D 



and 



Qk = f ■ 



for k < mo , 



(57) 



(58a) 



Q k = f + 1- , for k > m D , (58b) 

where /j, / and ran are set by the user (standard reference 
values are /j = 1, / = and mo = 0, 1). 

Besides being particularly dissipative, this method has the 
important advantage that it can be adapted to spherical sym- 
metry by modifying (52) as 



S D (^QdrF^dr^D. 



(59) 



In this way, the condition (49) is satisfied without the need of 
adding corrective terms. In EDGES, the term (59) is added us- 
ing an operator-splitting technique, so that its use corresponds 
to the application of a linear filter. This is discretized with 
two different techniques: locally on each element, or glob- 
ally, using an interior penalty technique. We call the latter ap- 
proach "interior-penalty spectral-viscosity" (IPSV) stabiliza- 
tion. While the first approach is completely local, the second 
one is able to relax also the jump terms between the elements. 
The main drawback of this method, however, is that we have 
found that the quality of the results can be quite sensitive to 
the tuning of [i and m D , which are necessarily problem de- 
pendent. 

The third method implemented in EDGES is usually re- 
ferred to as "spectral filtering" (see e.g. [55]). The idea behind 
this technique is to filter the numerical solution, or its time 
derivatives, with a low-pass filter in order to remove high- 
frequency components and keep the Gibbs oscillations under 
control. In the context of Legendre pseudospectral methods, 
this filtering, which we indicate with To, reads 



D 



[F D u] (x) = a ( Jj ) "fc L k 0) 



k=0 



(60) 



where Uk are the coefficients of the Legendre expansion of u 
and a(i]) is a filter function of order p in the Vandeven's sense, 
that is a function a E C P (R + ; [0, 1]) such that 

■*(0) = 1; (61) 
•<7 (fc) (0)=0, for k = 1,2,... ,p- 1; (62) 
• a(r]) = , for 77 > 1 ; (63) 

•<r (fc) (l)=0, for k = 1,2,... ,p- 1; (64) 

where f^ p \x) denotes the p-th derivative of / in x. These 
schemes were studied by Hesthaven and Kirby [56], who 
showed that if <r(-) is a filter of order p, and u e C p (— 1, 1), 
then 



\u{x) - T D u(x)\ < Z) 1 - p ||7t^|| L 2 ( _ 



1,1) 



(65) 



so that, if p <~ D, the filter does not spoil the spectral accuracy 
of the method. 

The effects of filtering on the stability of the numerical 
method can be intuitively understood from the fact that some 
filters, in particular the exponential filter of order p, 



17(77) = ex P 



(66) 



can be interpreted as equivalent to the use of numerical dif- 
fusion of order p and strength /j [57]. At present, however, a 
mathematical understanding of the impact of filtering on the 
accuracy and the stability of Legendre pseudospectral meth- 
ods is still lacking [56]. 

EDGES provides a number of different filter functions, but 
the most used are the exponential filter (66), which is only an 
approximate filter, but which is very popular due to its flexi- 
bility [55], and the "Erfc-Log filter" proposed by Boyd [58], 

, , 1 f L 1/2/1 1 1\ / - log[l - 4fo - lffl l 
a(7 7 ) = -crfc|2p/ (M-jjy 4(3-1/2)' /' 

(67) 

which is a semi-analytic approximation of the Vandeven's fil- 
ter, 



(7(7?) = 1 



r(2p) 
r(p)2 



/ [t(l-tT 
Jo 



dt. 



(68) 



where T is the Euler Gamma function. Note that it is custom- 
ary to use filters of order p <~ 2D or larger when expecting 
regular solutions and 2D + 1 is a reasonable value for a strong 
filter such as Erfc-Log. 

The main limitation of spectral filtering is that, being 
a high-order linear method, it cannot completely remove 
the Gibbs oscillations. This is a consequence of the well- 
known Godunov theorem stating that no linear monotonicity- 
preserving method exists of second-order or higher (see 
e-g- [6])- The idea is, instead, to allow for small oscillations 
of the numerical solution at the location of shocks, while pre- 
venting them from growing out of control. In our numerical 
experience we found that spectral filtering is a robust and vi- 
able alternative to slope limiting in the shock-tube case, where 
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Figure 2. Solution for the rest-mass density of a relativistic shock 
tube with different stabilization techniques. The solid black line rep- 
resents the exact, analytic solution, while the red crosses and the blue 
squares refer to the use of the two different DLT filters, i.e. the Erfc- 
Log and the exponential filter, respectively. The blue triangles and 
the violet circles refer instead to the use of the IPSV filter and of 
high-order slope limiters, respectively. To avoid excessive cluttering, 
only a data point every five is shown. The results were obtained on a 
grid of N = 200 elements with a polynomial of degree D = 3. 



its use often results in a sharper resolution of the discontinu- 
ities, and, at the same time, it is well-behaved in the spheri- 
cally symmetric case. Spectral filtering seems also to be much 
less sensitive than the spectral viscosity methods to the tuning 
of the parameters of the filter functions. Finally, we point out 
that, in contrast to IPSV and slope limiters, this technique is 
completely local in the sense that its application in a given 
element does not require the knowledge of the solution else- 
where. For these reasons, most of the results that we will show 
in the next Sections were obtained with the use of spectral fil- 
tering. 

Before doing that, however, and to illustrate the differ- 
ence between the various stabilization techniques, we show 
in Fig. 2 a comparison of the results obtained using different 
filters, while keeping all the other discretization parameters 
fixed. In particular, we show the results obtained in the case 
of a relativistic shock tube in flat spacetime. The details of 
the physical setup will be given in Sect. IV A, here instead we 
focus on the effects of the different stabilization techniques. 
As anticipated we find that slope-limiting, even in their high- 
order variant, result in a much larger smearing of the shock 
front, then all the other methods. The 9th order Erfc-Log fil- 
ter yields a much sharper resolution of the shock and a stable 
evolution, even tough, as discussed before, filtering does not 
have such a strong mathematical basis, while slope-limiting is 
known to yield a total-variation diminishing in mean (TVDM) 



scheme [19]. The best overall results are the ones obtained 
with the second-order exponential filter, for which we used a 
small strength factor, \i = 1, and the IPSV technique with 
/j = 1 and Qk — 1. We remark that the results obtained with 
these last two methods are very similar since, as discussed be- 
fore, the action of an exponential filter is roughly equivalent 
to that of a spectral viscosity. The results in Fig. 2 obtained 
on a grid of N — 200 elements with a polynomial of degree 
D = 3. 



D. Treatment of low-density regions 

The treatment of interfaces between vacuum region and 
fluid regions is one of the most challenging problems in Eu- 
lerian (relativistic) hydrodynamics codes (see e.g. [59-61]). 
The most commonly used approach to treat vacuum regions 
is to fill them with a low-density fluid, such that if a fluid 
element is evolved to have a rest-mass density below a cer- 
tain threshold, it is set to have a floor value and zero coordi- 
nate velocity [62, 63]. This approach works reasonably well 
and has been adopted by the vast majority of the relativistic - 
hydrodynamics codes. Nevertheless, the introduction of this 
"atmosphere" creates several numerical difficulties. First of 
all, this approach is gauge dependent and, for example, a star 
can accrete or lose mass from/to the atmosphere due to oscil- 
lations in the coordinate variables. Secondly, this approach 
often results in the introduction of errors at the surface of the 
star, which, on the one hand, usually have small influence on 
the overall dynamics because of the small momenta involved, 
but which, on the other hand are relatively large if compared 
with the magnitude of the involved quantities at the surface. 
Finally, and most importantly, it is possible that the algorithm 
governing the evolution of the floor can couple with the Gibbs 
oscillations, leading quickly to their amplification and destabi- 
lizing the scheme (After-all, the introduction of an atmosphere 
treatment is de-facto equivalent to the use of a boundary con- 
dition and this can very well lead to an unstable algorithm). 
For this reason this procedure could not be used in a straight- 
forward way in our code, but required a particular care. 

Other solutions, such as the use of the equations in La- 
grangian form, e.g. [64], the use of moving grids techniques 
e.g. [13], or the use of suitable limiters at the surface [59], do 
not suffer from these issues, but are restricted to the spheri- 
cally symmetric case. As our code is meant to be a prototype 
code to study the viability of DG methods for relativistic hy- 
drodynamics, the use of these techniques would have defeated 
our purpose. 

As a result, several different approaches were implemented 
and tested in EDGES to overcome the difficulties with the at- 
mosphere discussed above. The main idea behind all these ap- 
proaches was to use stronger stabilization techniques at the in- 
terfaces between fluid and vacuum regions. Unfortunately all 
these techniques performed quite poorly because they resulted 
in an unacceptable lowering of the resolution at the surface 
and thus in large numerical errors. Moreover, in the spheri- 
cally symmetric case, the situation is greatly worsened by the 
fact that these errors tend to build up coherently while travel- 
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Figure 3. Division of an element into control volumes for the 
spectral-volume method applied at the fluid-vacuum interface. The 
values of the solution u on the collocation points (filled squares) are 
interpreted as volume-averaged values of the solution in appropriate 
control volumes (hollow diamonds denote the boundaries). These 
values are then used to generate a reconstructed solution IZu using 
an adapted slope-limiter method. 



ing towards the center. For example, in the case of a neutron 
star, they tend to be amplified by a factor i?/r^l0 3 ,i?^10 
being the radius of the star and r ~ 0.01 being the location of 
the closest grid point to the center. 

Within EDGES a solution to this problem was eventually 
found with the use of a method that is, at the same time, 
oscillation-free and capable of obtaining high-enough resolu- 
tion at the surface. In particular, we derived a "sub-element 
method" approach that we later discovered to be very simi- 
lar to the spectral volumes (SV) method already suggested by 
Wang [65]. In this method, we first flag those elements in 
which the rest mass density falls below a certain threshold. 
Secondly we interpret the value of the solution on the collo- 
cation points of those elements as being the volume-averaged 
value of the solution on appropriated control volumes. Finally 
these values are evolved as in a FV scheme: a linear recon- 
struction with slope limiting is used to compute the value of 
the primitive variables at the interface of the control volumes 
and the HLLE approximate Riemann solver is used to com- 
pute the fluxes. 

The structure of the grid used for this spectral volume 
method is shown in Fig. 3. The control volumes are simply 
obtained by considering the points of the dual mesh as the 
cell interfaces, while their volume has to be corrected to take 
into account the weights of the Gaussian quadrature associ- 
ated with the primal grid. The linear reconstruction can be 
performed only for the interior control volumes, while for the 
control volumes at the boundary of each element we already 
have the value of the solution at one side of the control vol- 
ume, thus we are forced to reconstruct the solution there as a 
constant. In our tests we found that the "superbee" [3] limiter 
is the one which guarantees the best results among the ones 
that we tested. For this reason, all the results shown in the 
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Figure 4. First shock-tube test. Shown with solid lines are the ex- 
act solution for the rest-mass density (black line), the velocity (blue 
line), and the specific internal energy (red line). The correspond- 
ing numerical solution is represented by the black crosses, the blue 
triangles, and the red squares, respectively. To avoid excessive clut- 
tering we show only one point every 15 of the numerical solution. 
The results were obtained on a grid of N = 500 elements with a 
polynomial of degree D — 3. 



following make use of this limiter. 

The important advantage of this approach, with respect to a 
more traditional approach in which the troublesome elements 
are split into equal-size cells and a classical FV method is 
used, see e.g. [66], is that no interpolation is required in the 
switching from/to the discontinuous Galerkin method. The 
coupling between the two methods is also very natural and 
is done through the Riemann solvers between the elements. 
Again no special treatment is required to handle different type 
of elements: the DG ones and the SV ones. The main limita- 
tion of this approach is that, as we do not increase the number 
of degrees of freedom in the flagged elements, it has the effect 
of reducing the accuracy to second order where it is used. As 
it will be shown in Sect. IV D, this seems not to be limiting 
the accuracy of our code, probably because the use of the SV 
method is confined to small regions containing low-density 
fluid. In the case of a neutron star, for example, the use of the 
SV method is typically limited to one element containing the 
surface of the star. 

As a final remark we note that, as we are working with 
a nonuniform grid, the slope limiting method in its standard 
form, also used by us, is not guaranteed to be TVD and/or 
second order [67]. Again this seems not to be a problem in 
practice. 
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Figure 5. Rest-mass density profile for the second shock-tube test. 
Shown with different lines are the exact solution (solid black), the 
solution obtained with the ECHO code [69] (dotted blue). Both codes 
use 1000 cells/elements and algorithms at fifth order; EDGES also 
uses a exponential filter with fi = l,p = 2. 



Figure 6. The same as in Fig. 5, but with a magnification of the 
shock and a comparison of the results at different orders D of the 
polynomial representation of the solution. Note that the accuracy 
with which the shock location is founded increases with D 



IV. NUMERICAL TESTS 

In what follows we present the results obtained from some 
of the tests performed with EDGES. These tests have been 
chosen because they highlight the capabilities of our code in 
idealized settings, such as in the simulation of shock tubes 
and shock waves, and in more "astrophysically motivated" 
settings, such as in the simulation of the dynamics of spheri- 
cal stars and of the gravitational collapse to black hole. In all 
cases considered, the evolution has been made employing an 
ideal-fluid EOS to take into account for non-isentropic trans- 
formations, such as shock-heating effects. 



A. Shock tubes 

The first tests performed are shock-tube tests and, more 
specifically, we first present the results obtained in the case 
of two standard benchmarks for relativistic-hydrodynamics 
codes described in [68]. These are referred to as "blast- wave" 
problems 1 and 2 in [68], and are essentially one-dimensional 
flat-space Riemann problems. 

The first problem describes the propagation of a relativistic 
blast wave through a low-density medium. The solution ob- 
tained with EDGES, as well as the analytic solution, are shown 
in Fig. 4, where the solution is computed using polynomials 
of degree three (i.e. D = 3) and 500 elements (i.e. N — 500). 
The stabilization is obtained with the use of an Erfc-Log fil- 
ter of 9th order, which has been applied directly to the con- 
served variables in the regions where they fail to be monotone 



at every timestep. As it can be seen from the figure, the sta- 
bilization is strong enough to eliminate any oscillation, thus 
reducing the order of the method to the expected first-order 
near the discontinuity, yet with a very small smearing of the 
shock. 

The second shock-tube problem is similar to the first one, 
but is much more extreme. This benchmark is considered to 
be a challenging test also for modem HRSC codes [68]. The 
main difficulty is to handle the very high compression ratios, 
typical of relativistic hydrodynamics, produced in this case. 
In Fig. 5 we show a comparison of the results obtained with 
our code when using the exponential filter (66) at second- 
order (i.e. p — 2) with fi = 1.0, which is employed in the 
same way as the Erfc-Log filter in the previous test. Also 
shown for a comparison is the solution obtained with the high- 
order HRSC code ECHO [69], which was kindly provided 
by Olindo Zanotti. The ECHO code uses HLLE Riemann 
solver and monotonicity-preserving reconstruction at fifth or- 
der [70], while for EDGES we employed fifth-order polyno- 
mials 8 . Furthermore, the comparison is made when using the 
same number, 1000, of cells (for ECHO) and of elements (for 
EDGES); for both codes, in fact, we expect a convergence or- 
der (on smooth solutions) of the type Ax~ p , with Ax being 
the width of each cell/element and p being the order of the 
scheme. Note, however, that, while in ECHO high order is ob- 
tained with the use of data across multiple cells, in EDGES 



Note that the solution of the ECHO code used here does not make em- 
ploy the high-order DER flux reconstruction discussed in the Appendix of 
Ref. [69] and which would remove some of the oscillations. 
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this is attained with the use of a polynomial representation of 
the solution in each element. Thus the EDGES code actually 
uses around six times more grid points than the ECHO code. 
Nevertheless, we argue that this is the correct comparison be- 
cause, as mentioned before, the convergence properties of the 
two codes once the reconstruction procedure and the order of 
the polynomial representation is fixed, scale with the width 
of the cells/elements. This kind of comparison, where results 
obtained with discontinuous Galerkin and finite volume meth- 
ods are obtained using the same number of elements/cells, but 
different number of degrees of freedom has already been dis- 
cussed in the literature, see e.g. [71]. 

As can be seen from Fig. 5 the quality of the solution ob- 
tained by the two codes is comparable, with the solution ob- 
tained by the EDGES code being slightly closer to the exact 
position of the shock, but also showing signs of Gibbs os- 
cillations ahead of the shock. Overall, the exponential fil- 
ter can prevent the oscillations from growing and yields the 
best results when compared to the other flattening techniques 
discussed above. A more quantitative estimate of the qual- 
ity of the numerical solution can be obtained by looking at 
the ratio <7/<7 eX act, between the observed compression ratio 
and the analytic one [68]. We obtain a value of 0.96 with 
both ECHO and EDGES. At lower resolution, i.e. using only 
400 elements, we obtain, for EDGES, a value of 0.60 which is 
close to the ones reported for PPM at the same resolution: 
i.e. er/er cxac t — 0.54 — 0.69 [68]. We conclude therefore 
that in the case of discontinuous solutions, DG methods be- 
have similarly to FV methods with the same number of ele- 
ments/regions. 

It is worth stressing that while for ECHO the fifth order rep- 
resents the largest achievable, the EDGES code can perform 
this very severe test also with higher-order polynomials. This 
is shown in Fig. 6, where we report again the exact solution 
shown in Fig. 5 and compare it with the numerical one for 
D = 3,5,7 and D — 9. Of course the mathematical conver- 
gence in all of the different solutions is still only at first order, 
but it is remarkable to see that the code can handle such large 
shocks even when D = 9 and that the position of the shock 
and of the contact discontinuity become increasingly accurate 
as D is increased. Overall, the solution with D = 5 is the 
one yielding the closest match with the density right behind 
the shock and we have not considered orders higher than this 
in most of the subsequent tests. As a final remark we note 
that, to the best of our knowledge, this is the first time that a 
shock-tube test in relativistic hydrodynamics has been com- 
puted using a method of this order. 



B. Spherical shock reflection 

Another classical benchmark for relativistic - 
hydrodynamics codes is the spherical shock-reflection 
test. This is the relativistic version of the classical Noh test 
and its setup is described in detail in [68]. The initial data 
consist in a cold fluid converging at the center of the domain. 
As the fluid flows towards the center, a hot dense shell of 
matter is formed and a shock wave is propagated outwards. 
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Figure 7. The same as Fig. 4, but for the spherical shock-reflection 
test with a Lorentz factor of W = 1000. To avoid excessive clutter- 
ing in the main frame of the picture we plot only a point every five of 
the numerical solution, while the inset shows the solution on all the 
collocation points. 



In this situation it is possible to generate arbitrarily high com- 
pression ratios by increasing the Lorentz factor of the ingoing 
flow. This is a peculiarity of relativistic hydrodynamics and, 
again, the key point of the test is to assess the capability of the 
numerical methods to handle such strong density contrasts. 

This problem is also interesting to assess the quality of the 
correction procedure for the filtering in spherical symmetry. 
In particular, the quality of the results depends on the used 
correction procedure: they are of comparable quality when 
using the "dummy" and the "bubble correction" procedures 
(the former being slightly better), while the "intrinsic" correc- 
tion yields a filtering which is too weak to guarantee stability. 
This is not a surprise since with the intrinsic procedure we are 
not able to filter the coefficients of order less than two in the 
Legendre expansion, so that we cannot reduce the order of the 
method below the second even with a strong targeted filter. 

In Fig. 7 we present the results obtained for a flow with 
W = 1 000 solved using polynomials of degree three, on 
200 elements and with an Erf c -Log filter with p = 4 and a 
dummy correction. As shown in the figure, the filtering proce- 
dure is able to suppress any spurious oscillation and the shock 
front is captured within two elements, i.e. eight collocation 
points. The average relative error on the compression level 
is of 1.4%, comparable with the 2.2 % reported by Romero 
et al. [37] using FV with a minmod reconstruction on a grid 
with the same number of cells. Furthermore, as in Romero et 
al. [37], the error in e and p increases in the inner part of the 
solution. This is due to the well known wall-heating effect and 
is due to the small but nonzero dissipative and dispersive fea- 
tures of the numerical methods (see [72, 73] for an extended 
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Figure 8. L 2 -error in the spherical accretion test as a function of the 
number of elements N. Crosses and squares refer to the errors on Ci 
and C2, respectively, while the black, blue and red colors denote refer 
to the simulations using third, fourth, and fifth-order polynomials. 
Finally the dashed lines show the slopes associated with third, fourth, 
and fifth-order convergence. 



discussion of this effect). 

C. Spherical accretion onto a Schwarzschild black hole 

Having verified the capability of our code to handle shock 
waves, we next present some results concerning its accuracy 
in the case of regular solutions. In particular, we consider the 
case of spherical accretion onto a Schwarzschild black hole in 
the Cowling approximation (i.e. with a fixed spacetime). An 
analytic solution exists in the case of stationary flows and was 
first presented by Michel [74] and later used as a numerical 
test for other numerical codes, for example by [37]. This so- 
lution can be described in term of two constants: 

d=V=gr, (69) 

and 

c 2 = V^gT r t . (70) 

In our simulation we consider a spherical shell with ex- 
tent 3 — TV/40 < r < 20 in the spacetime metric of a 
Schwarzschild black hole of unit mass, while the initial condi- 
tions for the hydrodynamical variables describe a low-density 
fluid at rest with respect to an observer at infinity. At t = 
we start injecting higher-density fluid from the outer bound- 
ary and after a short transient the solution reaches stationarity, 
allowing us to measure C\ and C2 and to compare them with 
the analytic values fixed by the outer boundary condition. 



Figure 9. The same as in Fig. 8, but when the error is shown as a 
function of the polynomial degree when the number of elements N 
is kept fixed (black for N = 20 and red for N = 30). Note that an 
exponential convergence is measured. 



In Fig. 8 we show the L 2 -norm of the error for C\ and C2 
as a function of the number of elements N, and for different 
orders D of the polynomials used for the representation of the 
numerical solution over the single elements. The stabiliza- 
tion was obtained using an exponential DLT filter of order 2D 
and strength fj, = 1.0, corrected using the "dummy" strategy. 
Clearly, the errors computed from C\ and from C2 are basi- 
cally identical. More importantly the convergence order is the 
one that one would expect from the theory: third, fourth, and 
fifth-order using third, fourth, and fifth-order polynomials, re- 
spectively. 

A complementary measure of the convergence properties 
of the code is shown in Fig. 9, where we report the L 2 -norm 
of the error obtained after keeping fixed the number of ele- 
ments, N = 20, 30, but changing the order of the polynomi- 
als used for the representation of the solution. As expected, 
SDGM schemes behave in this case as multi-domain spectral 
methods [18], with the error decreasing exponentially with the 
polynomial degree. This is indeed what we observe in Fig. 9, 
where the errors decrease exponentially by almost three orders 
of magnitude for both C\ and C2, and at both resolutions. No 
sign of saturation appears in the error and we find it remark- 
able that the same method that was able to capture so sharply 
the discontinuous solutions in the previous tests, is also able 
to attain exponential convergence in smooth flows. 



D. Linear oscillations of spherical stars 

The results presented so far refer to situations that are some- 
how idealized and are meant mainly as a way to highlight the 
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Figure 10. Evolution of the central rest-mass density normalized to 
its initial value for a stable TOV and evolved with a polynomial of 
degree five. Different lines refer to different resolutions: solid black, 
dashed blue and dotted red for N = 75 , 100 and 125, respectively. 
Shown in the insets are the evolution at the initial time (i.e. t < 2 ms) 
and when the error at the surface is released triggering new small- 
scale oscillations (i.e. t ~ 45 ms). 



Figure 11. Power spectral density computed from the first 5000 Mq 
of the evolution of the central rest-mass density shown in Fig. 10. 
Different line types and colors refer to the different resolutions used 
(see legend). Shown as vertical dotted lines are the eigenfrequencies 
computed from linear perturbation theory. 



code's properties with respect to shock capturing, filtering cor- 
rection and accuracy on smooth solutions. We next present 
the results obtained in tests describing systems of more di- 
rect physical interest. Furthermore, in contrast to the previous 
tests, the matter will not be considered as a test-fluid and the 
spacetime will be properly evolved. 

A first interesting test is the study of linear oscillations of 
spherical stars. The initial data considered describes the equi- 
librium configuration of a self-gravitating fluid sphere, de- 
scribed by a polytropic equation of state, p = Kp v , and 
is obtained by integrating the Tolman-Oppenheimer-Volkoff 
(TOV) equations, as described by [75], and interpolating the 
result on the computational grid used by EDGES. The system 
is then evolved under the sole effects of numerical perturba- 
tions, mainly due to the interpolation errors in the generation 
of the initial data and to the interaction between the star and 
the atmosphere. Note that in order to avoid the presence of 
collocation points at r = 0, we stagger our numerical grid so 
that the r — point falls at the centre of the first element and 
we use only polynomials of odd degree. 

More specifically, we have considered a TOV constructed 
with a polytropic EOS having K — 100 and T = 2 and 
whose initial central density is p c — 1.28 x 10~ 3 . This 
could be taken to represent a stable, nonrotating, neutron 
star with gravitational mass M = 1.4 Mq and areal radius 
R = 9.6 Mq ~ 14.2 km. The degree of the polynomial basis 
is five and an Erfc-Log filter of order eleven is used to ensure 
a stable evolution. This filter is corrected with the "intrinsic" 



procedure outlined in Sect. Ill C 1 and is applied to the right- 
hand-side of the time stepping scheme. This results in a very 
weak stabilization algorithm whose main effect is to diffuse 
back into the atmosphere the numerical errors that would oth- 
erwise accumulate at the surface of the star and destabilize 
high-resolution runs. The latter, in fact, have very low intrin- 
sic viscosity and run for many more timesteps on timescales 
of 100 ms or longer (we recall that the dynamical timescale 
for this star is of the order of 0.7 ms). 

In Fig. 10 we show the evolution of the central density for 
three of these runs using different resolutions. They span a 
time frame of 20 000 Af Q ~ 100 ms, with the grid extending 
over the range < r < 15 and being composed by N = 
50, 75, 100, 125 and 150 elements. 

The dynamics is similar among the three runs, with the 
central density first exhibiting a "burst" with a variation of 
the order of 0.12 %, which is due to the fact that the star is 
"cut" by our interpolation algorithm of the initial data there 
where the density falls below the atmosphere threshold. This 
phase is magnified by the inset in the top of the figure. In 
a second stage the "burst" is quickly damped and the star 
starts to vibrate radially with its characteristic eigenfrequen- 
cies and eigenmodes under the effects of numerical perturba- 
tions, mainly at the surface of the star. In a third phase these 
oscillations are damped by the numerical viscosity, which 
depends on the resolution, and by the crude treatment used 
to represent the surface of the star. High-frequency modes 
are damped quickly while low-frequency modes have longer 
damping times. As a result, towards the end of the simulation 
we are left only with slowly-damped sinusoidal oscillations 
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associated with the F-mode. This is particularly evident in 
the N — 75 run, which has the largest numerical viscosity. 

The evolution is stable for all the resolutions that we have 
considered, but a careful examination of the behaviour of the 
central density reveals the existence of a fourth stage of the 
dynamics in some of the runs with lower resolution. More 
specifically, it is possible to note that during the third phase of 
the dynamics numerical errors accumulate at the surface of the 
star, producing small variations in the density profile near the 
atmosphere. When these variations are large enough so that 
they cannot be controlled by the employed flattening meth- 
ods, they are "released" and new energy is pumped into the 
high-frequency modes, starting a new phase in the dynamics. 
This happens at around 45 ms for the N = 75 run in Fig. 10, 
when a small high frequency component modulated by the F- 
mode suddenly appears, as shown by the inset at the bottom of 
the figure. The same happens around 80 ms for the N — 100 
run, while this phenomenon is not observed within the simula- 
tion time for higher resolution runs. Although the onset of this 
"energy release" at the surface can be delayed with the use of 
stronger filtering, it occurs only on a secular timescale and it 
does not affect the stability of the star. As a result, the central 
density keeps oscillating with a very small amplitude, of the 
order of 10 -3 %, even if, as discussed in Sect. IIIC2, pertur- 
bations at the surface tend to be greatly amplified when they 
reach the center as an artifact of the spherical symmetry. For 
this reason we believe that the energy release does not impact 
the quality of our simulations. 

A traditional measure of the accuracy of general-relativistic 
hydrodynamics codes is the comparison of the power spec- 
trum of the oscillations of TOVs against the values provided 
by linear perturbation theory. In Fig. 1 1 we make this com- 
parison by showing the power spectrum of the first 5 000 M Q 
of the evolution of the central density for different resolu- 
tions. The vertical dotted lines represent the eigenfrequencies 
computed from perturbation theory, which were kindly pro- 
vided us by Kentaro Takami and computed using the method 
described in [76]. We note that, even at the lowest resolu- 
tion, EDGES shows a perfect agreement between the observed 
proper frequencies and the perturbative ones for the F-mode 
and the first four overtones, HI, H2, iJ3 and HA. Further- 
more, as the resolution increases we are able to match more 
and more modes to the point that, with the N = 150 run, we 
match the first ten modes to a very good precision. We also 
note that the N = 150 run gives evidence of some nonlin- 
ear features in the spectrum, which shows considerable power 
also in its high-frequency part. 

A more quantitative measure of the convergence of the 
power spectrum is shown in Fig. 12, where we report as a 
function of the resolution the relative difference between the 
measured frequencies for the overtones H6, HI , H8 and H9, 
and the ones computed from perturbation theory. More specif- 
ically, the numerical frequencies were computed from the data 
during the first 5 000 Mq by using the procedure proposed 
in [77] . Namely, we computed the discrete Fourier transform 
(DFT) of the data with an Hanning window and used a three 
point interpolation of the power spectrum to correct for the 
incoherency error and determine the correct eigenfrequencies 
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Figure 12. Absolute relative difference between the estimated eigen- 
frequencies and the ones computed from linear perturbation theory, 
shown for different modes and as a function of the resolution. The 
lower-order modes are not considered because their error is smaller 
than the nominal one of the power spectral density and we report 
only the values for which a reliable measure of the error was possi- 
ble. Finally, indicated as black dashed lines are the slopes associated 
with fifth and sixth-order convergence. 



(this procedure is conceptually equivalent to a Lorentzian fit 
of the peaks of the DFT). Note that the lower-order modes 
are not shown because their precision is such that the error 
on those frequencies is well below the nominal uncertainty 
of the DFT and is basically dominated by the error on our 
measure. (The values of the error is not shown for those reso- 
lutions/frequencies that could not be determined in a reliable 
way: namely, the H8 and H9 modes for the N = 75 run.). 
Overall, we find that the measured convergence rate depends 
somehow on the specific mode, but it is compatible with a 
fifth-order convergence rate (cf. the fifth and sixth-order con- 
vergence rates which are shown as dashed lines). As we will 
comment later on, this results indicates that our largest source 
of error on the frequencies is not dominated by the low-order 
FV approach used at the surface of the star. Rather, the be- 
haviour in Fig. 12 shows that in the case of global quantities 
such as the oscillation frequencies of a TOV, the treatment of 
the surface is not the most critical element of the "error bud- 
get" in EDGES 9 . 

As a concluding consideration we note that in a recent work 
Cerda-Duran [78] has proposed to measure the numerical 



9 Note that the rapid decrease of the error in the estimate of H6 and H7 
at high resolution is most probably due to the nonlinear effects mentioned 
above which concentrate power in these higher-order modes, making them 
sharper and better resolved. 
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Figure 13. Estimated damping rates for various eigenmodes as a 
function of the resolution as computed by EDGES in ID and by 
Whisky in 3D. Black triangles, blue squares and red circles refer 
respectively to the F , HI and H2 modes, while filled/hollow points 
distinguish the estimates with EDGES from those with Whisky. Fi- 
nally the dashed lines show the estimated convergence order. 
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Figure 14. Evolution of the normalized central rest-mass density of 
an unstable TOV migrating towards the stable branch for a resolution 
of N = 150. The black solid line refers to a run employing an IPSV 
stabilization, while the red dashed line refers to a run employing a 
DLT filter; the blue horizontal line denotes the central rest-mass den- 
sity of the stable model associated with the initial configuration. The 
inset shows a magnification of the dynamics around 1 ms, at the first 
peak in the central density. 



bulk viscosity of general-relativistic hydrodynamical codes by 
looking at the damping time of the F-mode in the case of os- 
cillating TOVs. In particular, the rate of change with resolu- 
tion of the damping time (and which clearly increases with 
resolution) can be used as a measure of the convergence rate 
of the code. To explore this interesting suggestion we have 
computed the damping time by analysing more systematically 
the evolution of the central density. More specifically, using a 
sampling frequency of s» 1 M , we have built a discrete sig- 
nal which was then divided into chunks of 512 points with an 
overlap of 128 points. A DFT of each chunk was then com- 
puted using an Harming window and the power of the signal at 
the frequency associated with each mode was computed with 
a linear interpolation of the absolute value of the DFT. In this 
way we obtained an estimate of the energy in the mode for 
each time window. Finally we performed a least square fit of 
an exponential function to determine the damping time. (We 
have estimated the accuracy of this procedure to be better then 
10%, on the basis of tests performed with signals of known 
spectral properties.) 

In Fig. 13 we show the results of this measure when made 
on different modes. The data for the N — 50 , 75 and 
N = 100 runs have been obtained by windowing the evolu- 
tion in the interval < t < 5 000 to avoid spurious values due 
to the fourth phase of the dynamics described above, while for 
the N — 125 , 150 runs it was necessary to use the full data in 
the interval < t < 20 000 100 ms to obtain a reliable es- 
timate of the damping time for the F-mode. As a comparison, 
we also report the results obtained in 3D using the Whisky 



code [63, 79-81], when computed using a similar procedure. 
The simulations with Whisky were done using a PPM recon- 
struction and the HLLE Riemann solver for the hydrodynam- 
ics and using a fourth-order FD scheme for the evolution of 
the spacetime on a uniform grid also covering < r < 15 10 . 

We find that the convergence order of the damping time of 
the F-mode with Whisky is one, in agreement with the re- 
sults reported by [78] for the COCONUT code [26, 82]. In [78] 
it was argued that the reason why the order reduces to first 
is that the damping is active mainly at the surface of the star, 
where the numerical methods are only first order. The results 
found with EDGES show however a different behaviour, with 
the order of the damping being 3,2.5 and 2 for the F, HI and 
H 2 modes, respectively. Because the treatment of the surface 
in EDGES can be seen as a variant of the FV method used by 
both Whisky and COCONUT and is therefore only first-order 
accurate, our results suggest that the coefficient of the first- 
order (and surface-induced) error is much smaller than the co- 
efficient of the error coming from the high-order filtering pro- 
cedure, which then becomes the dominant source of damping. 
This is in agreement with the high convergence order found in 
the measurement of the eigenfrequencies (cf. Fig. 12). 

As a side comment we want to point out that we are able 



Note that, since Whisky and EDGES use different gauges, the two resolu- 
tions are only roughly equivalent. 
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Figure 15. Comparison in the evolution of the normalized central 
rest-mass density of the migrating star between a simulation in ID 
with EDGES (dot-dashed red line), and a simulation in 3D with 
Whisky (solid black line). Also in this case, the blue dashed line 
denotes the central rest-mass density of the stable model associated 
with the initial configuration and the inset shows a magnification of 
the second peak (see main text for details). 



Figure 16. Comparison of different stabilization techniques in the 
migration test. Shown in the two panels is the evolution of the central 
rest-mass density around the first peak. The left panel refers to a DLT 
filter on grids of N = 100 and N = 200 elements (blue solid and 
dashed red lines), while the right panel refers the IPSV stabilization 
technique with N = 100 and iV = 200 elements (blue solid and 
cyan long-dashed lines). Note that the DLT- stabilized run reaches 
convergences with a smaller number of elements. 



to attain higher-than-second order convergence in our results, 
because the time-step that we use is small enough so that the 
errors in the spatial discretization are dominating over the er- 
rors due to the time evolution, which is only second order. 



E. Nonlinear oscillations of spherical stars: the migration test 

As a direct extension of the analysis carried in the previous 
Section on linear oscillations, we next study large nonlinear 
oscillations which are produced as an equilibrium star model 
on the unstable branch of equilibria models migrates to the 
stable branch. This process has been used as a numerical test 
in 3D codes (see, e.g. [62, 63, 83, 84]), has been analyzed 
extensively in the past [38, 85, 86] and has gained special in- 
terest recently when it was shown that it exhibits a critical 
behaviour [87, 88]. 

Here, in particular, we have considered a TOV constructed 
with a poly tropic EOS having K — 100 and T = 2, central 
density p c — 7.0 x 10~ 3 , gravitational mass M = 1.49 M 
and areal radius R = 6 M@ ~ 8.8 km. The evolution, on the 
other hand, was made with an ideal-fluid EOS to properly take 
into account shock-heating effects. The numerical grid covers 
the region < r < 30 and the simulations reported used a 
polynomial representation of the solution of degree five. Two 
different stabilization techniques were used. A first one em- 
ployed an exponential filter of order six with // = 40, applied 
to the DLTs of the conserved variables and corrected with the 



"intrinsic" procedure. A second one used instead an IPSV sta- 
bilization with Qk = 1 — Sko and strength /i = 1. 

To trigger the migration on the stable branch, the star is 
perturbed with an outgoing velocity perturbation of the form 

v(x) = -\x 3 -3x\, x=^, (71) 

where A — 0.01. Under the effect this perturbation the star 
exhibits a violent expansion and migrates towards a new stable 
equilibrium configuration with a series of large-amplitude os- 
cillations. During these violent oscillations the exterior layers 
of the star tend to infall with higher velocity then the interior 
layers and this leads to the formation of shock waves that heat 
the neutron star matter and result in the ejection of a small 
portion of the material of the star. 

In Fig. 14 we show the evolution of the central density, nor- 
malized to its initial value, for two runs employing 150 ele- 
ments and different stabilization techniques. Because this test 
does not have an analytic counterpart, we have compared it the 
corresponding evolution performed with the Whisky code in 
a 3D simulation having N = 100 grid cells in each direction, 
PPM reconstruction and the HLLE Riemann solver. 

The comparison is offered in Fig. 15 and shows an ex- 
tremely good agreement. We note that because the two codes 
have intrinsically different initial truncation errors, the expan- 
sion phase in Whisky is slightly delayed with respect to the 
dynamics produced by EDGES, and we account for this dif- 
ference by shifting the data obtained by Whisky in order 
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Figure 17. Evolution of the lapse function at the center of an unstable 
TOV migrating towards the stable branch for two different runs using 
different stabilization techniques (black solid line for the IPSV sta- 
bilization, and red dashed line for the DLT filter) and N = 150. The 
inset shows a magnification of the dynamics around the first peak at 
1 ms and highlights that the solution is always smooth. 



to obtain an approximate alignment at the first minimum of 
the density. Notwithstanding the very good agreement, the 
main difference between the two solutions is the presence of 
high-frequency oscillations in the central density computed 
by EDGES, and in particular for the DLT runs and near the 
maximum value of the density. These "spikes" can be tracked 
down to the propagation of shock waves which are formed in 
the outer layers of the star during the collapse phase. They 
have initially a small compression factor, but they also tend to 
sum up coherently as they travel towards the center, thanks to 
the assumption of spherical symmetry, so that they result in 
strong variation of the density near the center of the star. This 
phenomenon has also been observed in other works in spheri- 
cal symmetry, employing standard finite-volume schemes, see 
e.g. [84], but is not usually observed in 2D or 3D simula- 
tions and, indeed, it is not present in the results obtained with 
Whisky. In particular Whisky employs Cartesian coordi- 
nates, which inevitably introduce preferred directions in the 
grid, thus preventing the spherical waves produced at the sur- 
face from interfering constructively and focussing towards the 
center. Another reason why these spikes are not observed in 
multi-dimensional solutions could be related to the larger nu- 
merical viscosity due to the necessity of using a much coarser 
resolution than in the ID case. This is confirmed by the fact 
that the use of the IPSV stabilization techniques, which allows 
us to introduce a much larger numerical dissipation, is able to 
greatly suppress this phenomenon. 

A more detailed comparison between the results obtained 
with the two stabilization techniques is shown in Fig. 16, 



where we offer a comparison between the evolution of the 
central density around the first bounce obtained with the two 
approaches and with different resolutions. Note that the DLT 
runs show signs of spikes in the density for all the resolutions 
(left panel), while the IPSV runs are much smoother (right 
panel). On the other hand, the results obtained with the DLT 
filters are already in a convergent regime and in fact the results 
obtained with N = 100 elements are very similar to those ob- 
tained by doubling the resolution. This is not the case for the 
solutions obtained with the IPSV filter, that seem to be only 
slowly approaching the DLT ones (cf. note in the right panel 
that the IPSV solutions approaches the DLT one as N goes 
from 100 to 200). The evidence that the simulations with the 
DLT filtering are already in a convergent regime and yet show 
spikes in the evolution can be interpreted, therefore, as a con- 
firmation that the latter are not a numerical effect, but rather 
an artifact of the symmetry, which leads to a focusing of the 
waves travelling towards the center. 

Fortunately, because these perturbations actually carry only 
a very small energy and are amplified by the symmetry, the 
evolution of the spacetime variables is totally unaffected. This 
is shown in Fig. 17, where we report the evolution of a partic- 
ularly representative metric quantity, namely, the lapse func- 
tion at the center. As can be seen from the figure, the lapse 
function shows no spikes or spurious oscillation and appears 
to be smooth with both the DLT and the IPSV filters, the main 
difference being the value of the minimum attained during the 
first bounce with the former stabilization technique. 



F. Gravitational collapse of unstable spherical stars 

As a final test we consider another classical testbed in 
general -relativistic hydrodynamics: namely, the gravitational 
collapse to black hole of an unstable TOV star. The problem 
of the gravitational collapse to a black hole has been already 
studied in great detail in a number of different conditions in- 
volving ID, 2D and 3D simulations, as well as different phys- 
ical conditions (see, e.g. [37, 38, 63, 85, 86, 89-95]), and has 
become a standard test of general-relativistic codes. For this 
reason, we will not discuss the details of the dynamics of the 
collapse and concentrate instead on the quality of the results 
obtained with EDGES. 

For the tests considered here we have evolved an unstable 
TOV built with a polytropic EOS with K = 100 and T = 2, 
having central density p c = 4.5 x 10~ 3 , gravitational mass 
M = 1.6 M & and areal radius R = 6.9 M ~ 10.2 km. 
In contrast to the migration test, the collapse is triggered by 
introducing an ingoing velocity perturbation of the type (71) 
with A = —0.01. The evolution of this system is studied on 
a grid covering < r < 15, staggered about the origin, using 
polynomials of degree five and employing an exponential fil- 
ter of order six, with intrinsic correction and strength fj, = 40, 
applied directly on the conserved variables at every time step. 

During the evolution and as expected, the central density 
of the star shows an exponential increase, which halts when 
the lapse function collapses to zero, signalling the formation 
of the black hole and "freezing" the evolution in the inner re- 
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Figure 18. Evolution of the lapse function at the center of an un- 
stable TOV collapsing to a black-hole for different resolutions on. 
The solid black line, the dotted blue line and the dashed red line are 
associated with runs employing a grid of 50, 75 and 100 elements, 
respectively. The inset shows a magnification of the dynamics around 
0.6 ms, highlighting some post-collapse dynamics in the lapse for the 
low-resolution run, which is however absent at higher resolutions. 



gions of the star. The maximum value attained by the central 
density does not have any physical meaning, partly because it 
depends entirely on the gauge conditions and partly because 
the gauge used will prevent its growth past a certain value. For 
this reason it is more interesting to look at the evolution of the 
lapse function at the center of the star, and which we show in 
Fig. 18. Note that for all the considered resolutions, the lapse 
function quickly collapses to zero and shows no appreciable 
subsequent evolution, with the exception of the lowest reso- 
lution case, for which a small growth appears again around 
0.6 ms. This sudden evolution is the result of errors coming 
from the surface of the frozen star, where the spatial slice is 
stretched due to the fact that we fixed our shift gauge condi- 
tion to zero and, thus, the metric functions present large spa- 
tial gradients. These errors, which cannot be compensated by 
simply increasing the resolution and are a shortcoming of the 
gauge used, propagate towards the inner regions of the frozen 
star and can induce a growth of the lapse, restarting the dy- 
namics. This is clearly a resolution-dependent effect which 
disappears quickly by increasing the resolution (see inset in 
Fig. 18). Apart from this late-time, and low -resolution dy- 
namics, the three curves are on top of each other, signaling 
a convergent regime, with a rate we measure to be slightly 
above the fourth order before the collapse of the lapse. After- 
wards the convergence order starts to slowly decay, due to the 
fact that, as the lapse collapses to zero, the equations are only 
weakly hyperbolic. At the end of the simulation the conver- 
gence order reaches a value which is > 1.5. 



In summary, the solution of the gravitational collapse to a 
black hole with EDGES has been straightforward and indeed 
a test which is much less demanding that the migration one, 
at least in terms of the hydrodynamical variables. The only 
difficulty encountered is a well known one and has to do with 
the fact that with such a gauge the collapse of the lapse is 
not compensated by any change in the shift vector, which is 
identically zero. As a result, as the black hole is produced it 
quickly produces a stretching of the coordinates at the location 
immediately outside the "frozen star". This, in turn, results in 
the development of strong gradients in the numerical solution 
which cannot be resolved indefinitely without suitable adap- 
tivity. 



V. CONCLUSIONS 

Numerical, relativistic hydrodynamics and MHD have seen 
a tremendous growth in the number and the quality of its re- 
sults after the introduction of modern high-resolution shock 
capturing schemes [3]. These schemes have been proven to 
be of central importance in the modelling of complex sys- 
tems involving strong gravity and/or high Lorentz factor. Yet, 
they suffer from important limitations that ultimately impact 
the accuracy of the obtained results [11]. For this reason the 
search for better numerical schemes is still on-going (and will 
always be). 

Discontinuous Galerkin schemes were developed to over- 
come the previously mentioned limitations of finite-volume 
and finite-difference schemes, while maintaining important 
properties of these schemes, such as conservation and shock 
capturing, that made them so successful in a number of appli- 
cations [19]. For this reason they are a natural candidate as 
an alternative to more traditional methods also in relativistic 
hydrodynamics. 

In this paper we have developed the necessary mathemat- 
ical framework needed for the application of discontinuous 
Galerkin schemes to relativistic hydrodynamics in curved 
spacetimes. In particular, we have presented both a mani- 
festly covariant weak formulation of relativistic hydrodynam- 
ics and a more traditional one obtained within a 3 + 1 split. We 
have then specialized the latter formulation to the spherically 
symmetric case and implemented it in a new one-dimensional 
relativistic hydrodynamical code, EDGES, which uses a high- 
order spectral discontinuous Galerkin method. 

The code was tested in a number of situations, including 
shock waves, spherical accretion, linear and nonlinear oscil- 
lations of relativistic spherical stars and the gravitational col- 
lapse of unstable stars to black holes. Our results show that 
discontinuous Galerkin methods are able to sharply resolve 
shock waves and, at the same time, attain very high, spectral, 
accuracy in the case of smooth solutions. For this reason they 
constitute an excellent alternative to classical finite-volume 
and finite-difference schemes for relativistic hydrodynamics, 
especially in those situations in which shock waves as well 
as small-scale features of the flows have to be resolved. In 
light of the promising prospects shown with these tests and 
of their affinity with a pseudospectral solution of the Einstein 
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equations, we anticipate that discontinuous Galerkin methods 
could represent a new paradigm for the accurate modelling in 
computational relativistic astrophysics. 

As a final remark we note that the success of these methods 
in multi-dimensional implementations will ultimately depend 
on the development of techniques such as local timestepping, 
/ip-adaptivity and load balancing [18], which are needed to 
take full advantage of the flexibility of these schemes [96]. 
This will represent the focus of our future work. 
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Figure 19. Schematic representation of a causal slice (J, S) contain- 
ing the element Q j . The future and past boundaries d ± , as well as 
the time-like or null-like boundaries d x flj of Slj are indicated with 
blue and red lines, respectively. Finally, the shaded region represent 
the causal slice. 



An example of such a slice of the triangulation is given by 
the shaded region in Fig. 19: a causal slice is basically a slice 
whose timelike spatial boundaries d + and d~ are parts of a 
Cauchy foliation of the spacetime. A causal slice (/, S) is 
said to be a minimal causal slice if it also satisfies the CFL- 
like condition 
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Appendix A: Causal triangulations 

In what follows we provide some basic definitions related 
with the causal structure of the triangulation. 

For every open set A, we introduce the notation d + A, to 
indicate the future boundary of A, that is the set of all the 
points, p € OA, such that 

j + ( P )nA = d), (Ai) 

J + {p) being the causal future of p. Analogously we call past 
boundary of A, d~ A, the set of all the points p e dA for 
which there exist an open neighborhood, U, such that 

J + (p)nUdA, (A2) 

where we have indicated with A the closure of A. Finally we 
will write d x A for dA \ [d+A U d~ A]. 

For a generic set E, instead, the above definitions are mod- 
ified by considering the interior of the set, E. In the case of an 
element of the grid, ilj, these definitions are useful to identify 
the regions of the boundary of flj where the characteristics are 
always ingoing, d~£lj, always outgoing, d + ilj, or of mixed 
nature, d x flj. These different parts of the causal slice (/, S) 
containing the element flj are shown in Fig. 19. 

We also define as slice of the triangulation, or simply slice, 
any tuple (I, S) where J C {1, 2, ... , N} and S = \J jeI flj 
is connected. We will say that a slice (/, S) is a causal slice if 

d x Scdn. (A3) 



that is, if the characteristics originating in each element, 
only intersect the boundaries of the element itself, 90, or the 
future boundaries of the other elements in the slice. 

Finally we say that a triangulation follows the causal struc- 
ture of the spacetime if it can be written as union of minimal 
causal slices. Intuitively a triangulation that follows the causal 
structure of the spacetime is one that can be sliced in minimal 
causal slices, which are grids associated with a Cauchy foli- 
ation of the spacetime and which satisfy a CFL condition on 
each hyper-surface of this foliation. 

If the triangulation {Qj}f = i follows the causal struc- 
ture of the spacetime, in the sense defined above, then we 
can introduce a family of successive slices {{Ik, Sk)}^ =0 , 
parametrized by a discrete time k, such that d~ Sk+i C 
d + Sk U dil and fl — [J k Sk- We can then construct a globally 
explicit, locally implicit, scheme by solving (9) on the suc- 
cessive slices. The reason is that, for all j e Ik+i, the fluxes 
across d~flj and d x flj are completely determined by the data 
on Sk and the boundary conditions, while the fluxes across 
d + Vlj do not couple f2 3 with other elements within Sk+i- 



Appendix B: The continuity equation in the D = 3 case 

We next present an example in which we show in the ex- 
plicit form the discontinuous Galerkin scheme for the simplest 
of the equations of relativistic hydrodynamics, namely the 
continuity equation, and with a polynomial of order D = 3. In 
particular we will show the explicit form of the semi -discrete 
equations for the collocated values in a given element Sj . In 
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the EDGES code, nodes, weights and interpolating polynomi- 
als are always computed numerically. However, for clarity, 
we are going to present, here, the analytic expressions for the 
considered case. 

The Legendre polynomial of degree three in [— 1, 1] is given 

by 

L 3(x) = -(5x 3 -3a;), 

so that the Gauss-Legendre-Lobatto quadrature points, the ze- 
ros of (1 — x 2 )dL 3 (x)/dx, are 

{zi}? = o = {-1, -l/\/5, 1/75,1}. 

The Lagrange polynomials associated with the quadrature 
points are 

l (x) = ^(1 - x)(-l + 5x 2 ) , 

h(x) = |(-1 + a;)(l + a;)(-l + s/bx) , 
o 

h{x) = -|(-1 + a;)(l + a;)(l + y/lx) , 
8 

l 3 (x) = ^(l + x)(-l + 5x 2 ), 
and the associated weights, lu, = J j U{x) Ax, 
W?=o = {1/6, 5/6, 5/6, 1/6}. 



Dropping now the element index j in the notation of 
(40), so that for any function /, the symbol /, will denote 
f(<Pj(xi)), and restricting ourselves to the continuity equa- 
tion in Sj we find 



tit 



= ^2w k rlX k v k V k 



fe=0 



w t \r 3 - r \ 



dh(x k ) 
da; 



(Bl) 



where J 7 ^ denotes the fluxes of V, between Sj and Sj±i, as 
obtained by the HLLE solver, i.e. 



^ _ A+ - A- 



and 



T+ - 

^ - A+ - A" 



+ A + A-(2? J -_i,3-X> J - ) 



A + (^+i,o^+i,o) - A (^-,3^,3) 



+ A + A-(D j+ i, -2?j,3) 



(B2) 



(B3) 



where vjj , Vjj are the values of v and V at the i-th colloca- 
tion point in the j-th element and A + and A - are the maximum 
and minimum characteristic speeds, respectively. 
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